load signac RDS object:

signac<-readRDS(paste0(projPath,"/H3K27ac_BinClusteredObject.RDS"))
# 
# metadata<-readRDS(paste0(projPath, "metadataObject.RDS"))

annotations<-readRDS("/scratch/kristian/scATAC/signac/EnsDb.Mmusculus.v79_getRangesAnnotation.RDS")

Read in peak lists from MACS2:

MACS2: pseudobulk broadpeak:

broadpeak<-read.table(paste0(projPath, "/macs2/BinClusters/broad/",ExptName,"_", ExptName,"_AllCellrangerFragments_peaks.broadPeak"), header = FALSE,stringsAsFactors=FALSE)
macs2Pseudobulk<-remove_blacklisted2Granges(peakTable=broadpeak)
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
signac$FrIP_macs2Pseudobulk<-FractionCountsInRegion(
  object = signac, 
  assay = 'bin5000',
  regions = macs2Pseudobulk
)

hist(signac$FrIP_macs2Pseudobulk, breaks=100)

VlnPlot(signac, features ="FrIP_macs2Pseudobulk" )

MACS2: Per BinCluster

files<-
  list.files(path=paste0(projPath, "/macs2/BinClusters/broad"),pattern =paste0(ExptName, "_BinCluster"), full.names = T)

files<-files[grep(pattern="_peaks.broadPeak",x=files)]

data <- files %>%
  # read in all the files, appending the path before the filename
  map(~ read.table( ., header = FALSE,stringsAsFactors=FALSE))

# binCluster_GR_white<-resize(binCluster_GR_white,width = 1000)
#############################
#### Collect all in one file:
#############################
binClustersDFlist<-lapply(data, function(x) remove_blacklisted2Granges(peakTable=x))
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
## 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
## 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
## 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
binClustersList<-GRangesList(binClustersDFlist)

macs2BinClusters<-unlist(reduce(binClustersList))
rtracklayer::export.bed(object = macs2BinClusters,paste0(projPath,"/macs2/BinClusters/broad/",ExptName,"_macs2BinClusters.bed"))


#plotting:

signac$FrIP_macs2BinClusters<-FractionCountsInRegion(
  object = signac, 
  assay = 'bin5000',
  regions = macs2BinClusters
)

hist(signac$FrIP_macs2BinClusters, breaks=100)

VlnPlot(signac, features ="FrIP_macs2BinClusters" )

Cellranger pseudobulk peaks:

cellrangerPeaks<-read.table(paste0(outsPath,"/peaks.bed"), header = FALSE,stringsAsFactors=FALSE)
cellrangerPseudobulk<-bed3_2white_Granges(peakTable=cellrangerPeaks)
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
signac$FrIP_cellrangerPseudobulk<-FractionCountsInRegion(
  object = signac, 
  assay = 'bin5000',
  regions = cellrangerPseudobulk
)

hist(signac$FrIP_cellrangerPseudobulk, breaks=100)

VlnPlot(signac, features ="FrIP_cellrangerPseudobulk" )

SEACR pseudobulk peaks:

seacrPeaks<-read.table(paste0(projPath,"/seacr/",ExptName,"_FilteredCells_SEACRpeaks_stringentTop0.01peaks.bed.stringent.bed"),
                       header = FALSE,stringsAsFactors=FALSE)

seacrPseudobulk<-bed3_2white_Granges(peakTable=cellrangerPeaks)
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
signac$FrIP_seacrPseudobulk<-FractionCountsInRegion(
  object = signac, 
  assay = 'bin5000',
  regions = seacrPseudobulk
)

hist(signac$FrIP_seacrPseudobulk, breaks=100)

VlnPlot(signac, features ="FrIP_seacrPseudobulk" )

SEACR: Per BinCluster

files<-
  list.files(path=paste0(projPath, "/seacr"),pattern =paste0( "BinCluster"), full.names = T)

files<-files[grep(pattern="SEACRpeaks_stringentTop0.01peaks.bed.stringent.bed",x=files)]

data <- files %>%
  # read in all the files, appending the path before the filename
  map(~ read.table( ., header = FALSE,stringsAsFactors=FALSE))

seacrBinClustersList<-lapply(data, function(x) bed6_ToWhite_GRanges(peakTable=x))
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
## 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
## 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
## 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
seacrBinClustersGRList<-GRangesList(seacrBinClustersList)

seacrBinClusters<-unlist(reduce(seacrBinClustersGRList))

rtracklayer::export.bed(object = seacrBinClusters,paste0(projPath,"/seacr/",ExptName,"_seacrBinClusters_AllInOne.bed"))


signac$FrIP_seacrBinClusters<-FractionCountsInRegion(
  object = signac, 
  assay = 'bin5000',
  regions = seacrBinClusters
)

hist(signac$FrIP_seacrBinClusters, breaks=100)

VlnPlot(signac, features ="FrIP_seacrBinClusters" )

ChipSeeker annotation of peaks

# binAnno <- annotatePeak(binClusters, TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene::TxDb.Mmusculus.UCSC.mm10.knownGene,level = "gene", , annoDb="org.Mm.eg.db")
# 
# vennpie(binAnno)
# plotAnnoPie(binAnno)
# upsetplot(binAnno)


#list of anno
filelist<- list(
  "cellrangerPseudobulk"=cellrangerPseudobulk, 
                "macs2Pseudobulk"=macs2Pseudobulk,
                "seacrPseudobulk"=seacrPseudobulk,
                "macs2BinClusters"=macs2BinClusters,
                "seacrBinClusters"=seacrBinClusters
                )

peakAnnoList <- lapply(filelist, annotatePeak, TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene::TxDb.Mmusculus.UCSC.mm10.knownGene,
                       tssRegion=c(-3000, 3000), verbose=FALSE)
names(peakAnnoList)<-names(filelist)

plotDistToTSS(peakAnnoList)

# genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
# names(genes) = sub("_", "\n", names(genes))
# compKEGG <- clusterProfiler::compareCluster(geneClusters   = genes,
#                          fun           = "enrichKEGG",
#                          pvalueCutoff  = 0.05,
#                          pAdjustMethod = "BH")
# dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")

genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)

### Choose which peak finding method for downstream analysis:

binClusters<-macs2BinClusters
#binClusters<-seacrBinClusters

signac$FrIP_binClusters<-FractionCountsInRegion(
  object = signac, 
  assay = 'bin5000',
  regions = binClusters
)

hist(signac$FrIP_binClusters, breaks=100)

VlnPlot(signac, features ="FrIP_binClusters" )

Create signac object from peak list

fragments <- CreateFragmentObject(
  path = paste0(outsPath, "/fragments.tsv.gz"),
cells=rownames(signac@meta.data),
  validate.fragments = FALSE
)
## Computing hash
#add new counts based on MACS2 peaks:
counts <- FeatureMatrix(
  fragments =fragments,
  features = binClusters,
  process_n = 50,
cells=rownames(signac@meta.data)
 )
## Extracting reads overlapping genomic regions
dim(counts)
## [1] 68396   636
signac[["peaks"]] <- CreateChromatinAssay(counts = counts,
                                          fragments = fragments@path,
                                          annotation = annotations
                                            )
## Warning in CreateChromatinAssay(counts = counts, fragments = fragments@path, :
## Overlapping ranges supplied. Ranges should be non-overlapping.
## Warning: Non-unique features (rownames) present in the input matrix, making
## unique
## Computing hash
DefaultAssay(signac) <- 'peaks'

signac
## An object of class Seurat 
## 454368 features across 636 samples within 3 assays 
## Active assay: peaks (68396 features, 0 variable features)
##  2 other assays present: bin5000, RNA
##  2 dimensional reductions calculated: lsi_, lsi_umap
# add the gene information to the object

### Save RDS object:
saveRDS(signac, "H3K27ac_peakMatrixObject_RAW.RDS")


# signac<-
#   subset(
#   x = signac,subset =   passed_filters<75000 &
#     FrIP_binClusters > 0.35 &
#     promoter_ratio >0.05
#   )

Binarize

signac[["binaryPeaks"]]<-signac[["peaks"]]
## Warning: Cannot add objects with duplicate keys (offending key: peaks_), setting
## key to 'binarypeaks_'
signac[["binaryBins"]]<-signac[["bin5000"]]
## Warning: Cannot add objects with duplicate keys (offending key: bin5000_),
## setting key to 'binarybins_'
signac<- BinarizeCounts(signac, assay=c("binaryPeaks", "binaryBins"))

Dimensionality reduction

Find top features from bins-> TF-IDF, SVD, UMAP

DefaultAssay(signac) <- 'peaks'

signac <- FindTopFeatures(signac, min.cutoff = "q5",verbose = T)
signac
## An object of class Seurat 
## 886758 features across 636 samples within 5 assays 
## Active assay: peaks (68396 features, 68396 variable features)
##  4 other assays present: bin5000, RNA, binaryPeaks, binaryBins
##  2 dimensional reductions calculated: lsi_, lsi_umap
signac <- RunTFIDF(signac)
## Performing TF-IDF normalization
## Warning in RunTFIDF.default(object = GetAssayData(object = object, slot =
## "counts"), : Some features contain 0 total counts
signac <- RunSVD(signac,
               reduction.key = 'LSI_',
               reduction.name = 'lsi_')
## Running SVD
## Scaling cell embeddings
DepthCor(signac, reduction = "lsi_")

RunUMAP

Excluding dimension #1 since it correlates with seq depth.

signac <- RunUMAP(signac, reduction = 'lsi_', dims = 2:ndim, reduction.name = 'lsi_umap',n.neighbors=15)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 16:08:01 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:08:01 Read 636 rows and found 29 numeric columns
## 16:08:01 Using Annoy for neighbor search, n_neighbors = 15
## 16:08:01 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:08:01 Writing NN index file to temp file /tmp/RtmpNLnAOc/filec5ee54cd2f3
## 16:08:01 Searching Annoy index using 30 threads, search_k = 1500
## 16:08:01 Annoy recall = 100%
## 16:08:04 Commencing smooth kNN distance calibration using 30 threads
## 16:08:08 Initializing from normalized Laplacian + noise
## 16:08:08 Commencing optimization for 500 epochs, with 11512 positive edges
## 16:08:12 Optimization finished
#signac <- RunUMAP(signac, reduction = 'lsi_', dims = 2:30, reduction.name = 'lsi_umap', n.start=1000,spread=0.24, min_dist=1, n.neighbors = 15)


p_bins<-DimPlot(signac, reduction = 'lsi_umap', label =T, label.box = T, repel = T, group.by = paste0("bin5000_snn_res.",BinRes))+
  ggtitle("Clusters from bin5000")+
  NoLegend()
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.
p_bins

FeaturePlot(signac, features="enhancer_region_ratio" ,min.cutoff = "q5",split.by = "is_cellranger_cell_barcode")

FeaturePlot(signac, features=c("nCount_peaks"),reduction = "lsi_umap" )

FeaturePlot(signac, features=c("FrIP_binClusters"),reduction = "lsi_umap" )

FeaturePlot(signac, signac@assays$peaks@var.features[1:4], slot="counts")

DefaultAssay(signac) <- 'peaks'

signac <- FindNeighbors(object = signac, reduction = 'lsi_', dims =2:ndim, force.recalc = T)
## Computing nearest neighbor graph
## Computing SNN
signac <- FindClusters(object = signac, verbose = FALSE, algorithm = 3, resolution=PeakRes)
## Warning: UNRELIABLE VALUE: Future ('future_lapply-1') unexpectedly generated
## random numbers without specifying argument '[future.]seed'. There is a risk that
## those random numbers are not statistically sound and the overall results might
## be invalid. To fix this, specify argument '[future.]seed', e.g. 'seed=TRUE'.
## This ensures that proper, parallel-safe random numbers are produced via the
## L'Ecuyer-CMRG method. To disable this check, use [future].seed=NULL, or set
## option 'future.rng.onMisuse' to "ignore".
 # table("binclusters"=paste0("signac@meta.data$bin5000_snn_res.",BinRes),"peak_cluster"=signac$seurat_clusters)
 # 

p_clusters<-DimPlot(signac, reduction = 'lsi_umap', label =T, label.box = T, repel = T, group.by = "seurat_clusters") +ggtitle("Clusters based on peaks")+NoLegend()
p_clusters

# p_dimplotSamples<-DimPlot(signac, reduction = 'lsi_umap', group.by="hash.ID", label =F, label.box = F, repel = T) +ggtitle("Samples based on peaks")
# p_dimplotSamples


 VlnPlot(signac, features="promoter_ratio", split.by="is_cellranger_cell_barcode")+ 
   VlnPlot(signac, features="logUMI", split.by="is_cellranger_cell_barcode")+
   plot_layout(guides = "collect")& theme(legend.position = 'bottom')
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

p_bins+p_clusters+plot_annotation(title="UMAP from peak matrix",theme =  theme(plot.title = element_text(hjust = 0.5)))

Make new bam-file for each cluster

cells<-
  signac@meta.data %>% 
  #as_tibble() %>%
  rownames_to_column("barcode") %>%
    filter(is__cell_barcode=="1") %>% 
  dplyr::select(barcode, seurat_clusters)

cells$seurat_clusters<-paste0("PeakCluster",cells$seurat_clusters)

write_tsv(x = cells, file=paste0(projPath,"/bam/FilteredCellBarcodes_PeakClusters.tsv"), col_names = FALSE)
VlnPlot(signac,
            features = c("nCount_bin5000",
                         "nFeature_bin5000",
                         "nCount_HTO",
                         "passed_filters",
                        #"mitochondrial",
                        "duplicate",
                         "enhancerPercentage",
                         "nCount_RNA",
                         "nFeature_RNA", 
                      #   "HTO_margin",
                         "blacklist_region_fragments","FrIP_binClusters"
                         )
)
## Warning in FetchData(object = object, vars = features, slot = slot): The
## following requested variables were not found: enhancerPercentage

  table(signac$seurat_clusters, signac$is_cellranger_cell_barcode)
##    
##       0   1
##   0  73 210
##   1   1 186
##   2   3 163
markers<-
  c(
#"Tbr1", 
"Slc17a6","Satb2","Nr5a1","Fezf1",
 "Slc32a1","Otp","Dlk1","Adcyap1",
"Agrp", "Pomc", "Npy", "Sst",
"Lepr",
"Meis2",
"Zfhx4", 
"Ghrh", 
#"Oxt",
#"Cck", 
"Ebf3",
#"Gfap",
"Nfib", "Esr1"
#"Rfx4", "Mog", "Alx4", "Gmds"
  )


DefaultAssay(signac) <- 'RNA'

p_markers<-FeaturePlot(
  object = signac,
  features = c(markers, "FrIP_binClusters"),
  pt.size = 0.1,
  max.cutoff = 'q90',
  ncol = round(length(markers)/5)
) & NoLegend() & NoAxes()

p_markers

p_vln_markers<-VlnPlot(signac, features=markers, group.by="seurat_clusters")
p_vln_markers

Save RDS object:

saveRDS(signac,paste0(projPath,  "/H3K27ac_PeaksClusteredObject.RDS"))

Find DA peaks

DefaultAssay(signac) <- 'peaks'
Idents(signac) <- "seurat_clusters"
da_peaks <- FindAllMarkers(
  object = signac,
  min.diff.pct = 0.1, #0.005
  logfc.threshold = 0.1,
  test.use = 'LR'
)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
# 
# da_peaks12<-FindMarkers(signac, ident.1 = "1", ident.2 = "3",
#   min.diff.pct = 0.1,
#   logfc.threshold = 0.1,
#   test.use = 'LR'
# )                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 

da_peaksGR <- makeGRangesFromDataFrame(separate(da_peaks, col=gene, sep="-",  c("Chr", "Start", "End")))
elementMetadata(da_peaksGR)<-da_peaks
names(elementMetadata(da_peaksGR))[7]<-"peak"

blacklist.df = read_tsv(file="/projects/jph712/CUTnTag/H3K27ac_ARC_scCUTnTAG_B/deepTools/mm10-blacklist.v3.bed" ,col_names = F) 
## 
## -- Column specification ----------------------------------------------------------------------------------------------------------------------------------------------
## cols(
##   X1 = col_character(),
##   X2 = col_double(),
##   X3 = col_double()
## )
blacklist.gr=GenomicRanges::GRanges(seqnames=blacklist.df$X1,
               ranges=IRanges(start=as.numeric(blacklist.df$X2),
                              end=as.numeric(blacklist.df$X3)))
simpleGR_diff<-GenomicRanges::setdiff(da_peaksGR, blacklist.gr, ignore.strand=TRUE)
da_peaksGR_white<-subsetByOverlaps(da_peaksGR,simpleGR_diff, ignore.strand=TRUE)

filtered_regions.df<- 
  da_peaksGR_white@elementMetadata %>% 
  as_tibble() %>% 
  filter(p_val_adj<0.1 & avg_logFC>0.25 |p_val_adj<0.1 &  -0.25>avg_logFC) %>% 
  distinct(peak, .keep_all=T) %>%
  group_by(cluster) %>% 
  arrange(avg_logFC) %>%
  dplyr::select(peak, cluster, avg_logFC,p_val_adj, everything())
DefaultAssay(signac) <- 'peaks'

####top in each cluster
closest <- ClosestFeature(signac,
  regions = filtered_regions.df$peak,
  )

filtered_regions.df<-left_join(filtered_regions.df, closest, by=c("peak"= "query_region"))
filtered_regions.df<-left_join(filtered_regions.df,
as_tibble(binClusters) %>%
 unite("peak", c(seqnames, start, end), sep = "-", remove = TRUE)
                               )
## Joining, by = "peak"
library(writexl)
writexl::write_xlsx(filtered_regions.df, path = "DiffAcetylatedPeaks_FindAll.xlsx", col_names = TRUE)
########

Plot cov plots:

covplot<-CoveragePlot(signac, 
             filtered_regions.df %>%
               group_by(cluster) %>%
               slice_min(p_val_adj, n=3) %>%
               pull(peak),
             extend.upstream = 5000, extend.downstream = 5000,
             group.by ="seurat_clusters",peaks = T,assay = "peaks" # not enough - needs to be default assay.
             ) 

suppressWarnings(print(covplot))

TilePlot(
  object = signac,
  region = "chr2-58289132-58290370", #Acv1r 
  tile.cells = 100,
  group.by = "seurat_clusters",order.by = "total",
  extend.upstream = 5000, extend.downstream = 5000
  )

goi<- c("chr15-66017872-66018112",
        "chr7-141678663-141678913",
        "chr3-37350491-37350692",
        "chr2-143343751-143344155",
        "chr13-51903198-51903570",
        "chr13-94856782-94881877"
)


CoveragePlot(signac, region = goi,
  group.by = "seurat_clusters", 
  extend.upstream = 5000, extend.downstream = 5000, ranges = binClusters, peaks=F)#,tile = T, tile.cells = 100)
## Warning: Removed 2 rows containing missing values (position_stack).
## Warning: Removed 12 rows containing missing values (geom_segment).
## Warning: Removed 12 rows containing missing values (geom_rect).
## Warning: Removed 18 rows containing missing values (geom_segment).
## Warning: Removed 18 rows containing missing values (geom_rect).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_rect).
## Warning: Removed 55 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 2 rows containing missing values (position_stack).
## Warning: Removed 10 rows containing missing values (geom_segment).
## Warning: Removed 10 rows containing missing values (geom_rect).
## Warning: Removed 19 rows containing missing values (geom_segment).
## Warning: Removed 19 rows containing missing values (geom_rect).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_rect).
## Warning: Removed 44 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (position_stack).

CoveragePlot(signac, 
               region = "chr3-119195567-119196174",
  group.by = "seurat_clusters", 
  extend.upstream = 5000, extend.downstream = 5000)
## Warning: Removed 21 rows containing missing values (geom_segment).
## Warning: Removed 21 rows containing missing values (geom_rect).
## Warning: Removed 38 rows containing missing values (geom_segment).
## Warning: Removed 38 rows containing missing values (geom_rect).
## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_rect).
## Warning: Removed 43 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_text).

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.3.so
## 
## locale:
##  [1] LC_CTYPE=C                 LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] writexl_1.3.1                            
##  [2] future_1.20.1                            
##  [3] org.Mm.eg.db_3.12.0                      
##  [4] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
##  [5] GenomicFeatures_1.42.1                   
##  [6] AnnotationDbi_1.52.0                     
##  [7] Biobase_2.50.0                           
##  [8] ChIPseeker_1.26.0                        
##  [9] GenomicRanges_1.42.0                     
## [10] patchwork_1.1.0                          
## [11] GenomeInfoDb_1.26.2                      
## [12] IRanges_2.24.1                           
## [13] S4Vectors_0.28.1                         
## [14] BiocGenerics_0.36.0                      
## [15] Signac_1.1.0                             
## [16] Seurat_3.2.2                             
## [17] forcats_0.5.0                            
## [18] stringr_1.4.0                            
## [19] dplyr_1.0.2                              
## [20] purrr_0.3.4                              
## [21] readr_1.4.0                              
## [22] tidyr_1.1.2                              
## [23] tibble_3.0.4                             
## [24] ggplot2_3.3.2                            
## [25] tidyverse_1.3.0                          
## [26] rmarkdown_2.5                            
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.1                         
##   [2] SnowballC_0.7.0                        
##   [3] rtracklayer_1.50.0                     
##   [4] GGally_2.0.0                           
##   [5] bit64_4.0.5                            
##   [6] knitr_1.30                             
##   [7] irlba_2.3.3                            
##   [8] DelayedArray_0.16.0                    
##   [9] data.table_1.13.4                      
##  [10] rpart_4.1-15                           
##  [11] RCurl_1.98-1.2                         
##  [12] AnnotationFilter_1.14.0                
##  [13] generics_0.1.0                         
##  [14] cowplot_1.1.0                          
##  [15] RSQLite_2.2.1                          
##  [16] shadowtext_0.0.7                       
##  [17] RANN_2.6.1                             
##  [18] bit_4.0.4                              
##  [19] enrichplot_1.10.1                      
##  [20] spatstat.data_1.5-2                    
##  [21] xml2_1.3.2                             
##  [22] lubridate_1.7.9.2                      
##  [23] httpuv_1.5.4                           
##  [24] SummarizedExperiment_1.20.0            
##  [25] assertthat_0.2.1                       
##  [26] viridis_0.5.1                          
##  [27] xfun_0.19                              
##  [28] hms_0.5.3                              
##  [29] evaluate_0.14                          
##  [30] promises_1.1.1                         
##  [31] fansi_0.4.1                            
##  [32] progress_1.2.2                         
##  [33] caTools_1.18.0                         
##  [34] dbplyr_2.0.0                           
##  [35] readxl_1.3.1                           
##  [36] igraph_1.2.6                           
##  [37] DBI_1.1.0                              
##  [38] htmlwidgets_1.5.2                      
##  [39] reshape_0.8.8                          
##  [40] ellipsis_0.3.1                         
##  [41] RSpectra_0.16-0                        
##  [42] backports_1.2.0                        
##  [43] biomaRt_2.46.0                         
##  [44] deldir_0.2-3                           
##  [45] MatrixGenerics_1.2.0                   
##  [46] vctrs_0.3.5                            
##  [47] ensembldb_2.14.0                       
##  [48] ROCR_1.0-11                            
##  [49] abind_1.4-5                            
##  [50] withr_2.3.0                            
##  [51] ggforce_0.3.2                          
##  [52] BSgenome_1.58.0                        
##  [53] checkmate_2.0.0                        
##  [54] sctransform_0.3.1                      
##  [55] GenomicAlignments_1.26.0               
##  [56] prettyunits_1.1.1                      
##  [57] goftest_1.2-2                          
##  [58] DOSE_3.16.0                            
##  [59] cluster_2.1.0                          
##  [60] lazyeval_0.2.2                         
##  [61] crayon_1.3.4                           
##  [62] labeling_0.4.2                         
##  [63] pkgconfig_2.0.3                        
##  [64] tweenr_1.0.1                           
##  [65] nlme_3.1-150                           
##  [66] ProtGenerics_1.22.0                    
##  [67] nnet_7.3-14                            
##  [68] rlang_0.4.9                            
##  [69] globals_0.14.0                         
##  [70] lifecycle_0.2.0                        
##  [71] miniUI_0.1.1.1                         
##  [72] BiocFileCache_1.14.0                   
##  [73] modelr_0.1.8                           
##  [74] rsvd_1.0.3                             
##  [75] dichromat_2.0-0                        
##  [76] cellranger_1.1.0                       
##  [77] polyclip_1.10-0                        
##  [78] matrixStats_0.57.0                     
##  [79] lmtest_0.9-38                          
##  [80] graph_1.68.0                           
##  [81] Matrix_1.2-18                          
##  [82] ggseqlogo_0.1                          
##  [83] boot_1.3-25                            
##  [84] zoo_1.8-8                              
##  [85] reprex_0.3.0                           
##  [86] base64enc_0.1-3                        
##  [87] ggridges_0.5.2                         
##  [88] png_0.1-7                              
##  [89] viridisLite_0.3.0                      
##  [90] bitops_1.0-6                           
##  [91] KernSmooth_2.23-18                     
##  [92] Biostrings_2.58.0                      
##  [93] blob_1.2.1                             
##  [94] qvalue_2.22.0                          
##  [95] parallelly_1.21.0                      
##  [96] jpeg_0.1-8.1                           
##  [97] scales_1.1.1                           
##  [98] memoise_1.1.0                          
##  [99] magrittr_2.0.1                         
## [100] plyr_1.8.6                             
## [101] ica_1.0-2                              
## [102] gplots_3.1.1                           
## [103] zlibbioc_1.36.0                        
## [104] scatterpie_0.1.5                       
## [105] compiler_4.0.2                         
## [106] RColorBrewer_1.1-2                     
## [107] plotrix_3.7-8                          
## [108] fitdistrplus_1.1-3                     
## [109] Rsamtools_2.6.0                        
## [110] cli_2.2.0                              
## [111] XVector_0.30.0                         
## [112] listenv_0.8.0                          
## [113] pbapply_1.4-3                          
## [114] ps_1.5.0                               
## [115] htmlTable_2.1.0                        
## [116] Formula_1.2-4                          
## [117] MASS_7.3-53                            
## [118] mgcv_1.8-33                            
## [119] tidyselect_1.1.0                       
## [120] stringi_1.5.3                          
## [121] yaml_2.2.1                             
## [122] GOSemSim_2.16.1                        
## [123] askpass_1.1                            
## [124] latticeExtra_0.6-29                    
## [125] ggrepel_0.8.2                          
## [126] grid_4.0.2                             
## [127] VariantAnnotation_1.36.0               
## [128] fastmatch_1.1-0                        
## [129] tools_4.0.2                            
## [130] future.apply_1.6.0                     
## [131] rstudioapi_0.13                        
## [132] foreign_0.8-80                         
## [133] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [134] lsa_0.73.2                             
## [135] gridExtra_2.3                          
## [136] farver_2.0.3                           
## [137] Rtsne_0.15                             
## [138] ggraph_2.0.4                           
## [139] rvcheck_0.1.8                          
## [140] digest_0.6.27                          
## [141] BiocManager_1.30.10                    
## [142] shiny_1.5.0                            
## [143] Rcpp_1.0.5                             
## [144] broom_0.7.2                            
## [145] later_1.1.0.1                          
## [146] RcppAnnoy_0.0.17                       
## [147] OrganismDbi_1.32.0                     
## [148] httr_1.4.2                             
## [149] ggbio_1.38.0                           
## [150] biovizBase_1.38.0                      
## [151] colorspace_2.0-0                       
## [152] rvest_0.3.6                            
## [153] XML_3.99-0.5                           
## [154] fs_1.5.0                               
## [155] tensor_1.5                             
## [156] reticulate_1.18                        
## [157] splines_4.0.2                          
## [158] uwot_0.1.9                             
## [159] RBGL_1.66.0                            
## [160] RcppRoll_0.3.0                         
## [161] spatstat.utils_1.17-0                  
## [162] graphlayouts_0.7.1                     
## [163] plotly_4.9.2.1                         
## [164] xtable_1.8-4                           
## [165] jsonlite_1.7.2                         
## [166] tidygraph_1.2.0                        
## [167] spatstat_1.64-1                        
## [168] R6_2.5.0                               
## [169] Hmisc_4.4-2                            
## [170] pillar_1.4.7                           
## [171] htmltools_0.5.0                        
## [172] mime_0.9                               
## [173] glue_1.4.2                             
## [174] fastmap_1.0.1                          
## [175] BiocParallel_1.24.1                    
## [176] codetools_0.2-18                       
## [177] fgsea_1.16.0                           
## [178] lattice_0.20-41                        
## [179] curl_4.3                               
## [180] leiden_0.3.6                           
## [181] gtools_3.8.2                           
## [182] GO.db_3.12.1                           
## [183] openssl_1.4.3                          
## [184] survival_3.2-7                         
## [185] munsell_0.5.0                          
## [186] DO.db_2.9                              
## [187] GenomeInfoDbData_1.2.4                 
## [188] haven_2.3.1                            
## [189] reshape2_1.4.4                         
## [190] gtable_0.3.0